Virtual-crystal approximation that works: Locating a composition phase boundary in 

Pb(Zri_,Ti,)03 
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We present a new method for modeling disordered solid solutions, based on the virtual crystal approximation 
(VGA). The VGA is a tractable way of studying configurationally disordered systems; traditionally, the po- 
tentials which represent atoms of two or more elements are averaged into a composite atomic potential. We 
have overcome significant shortcomings of the standard VGA by developing a potential which yields averaged 
atomic properties. We perform the VGA on a ferroelectric oxide, determining the energy differences between 
the high-temperature rhombohedral, low-temperature rhombohedral and tetragonal phases of Pb(Zri_a;Tia;)03 
at X = 0.5 and comparing these results to superlattice calculations and experiment. We then use our new method 
to determine the preferred structural phase at a; = 0.4. We find that the low-temperature rhombohedral phase 
becomes the ground state at a; = 0.4, in agreement with experimental findings. 
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Throughout condensed matter science, solid solutions pro- 
vide unique properties which are inaccessible in pure materi- 
als. Semiconductors and ferroelectric materials are two im- 
portant areas in which small changes in atomic composition 
dramatically change material properties. Theoretical exami- 
nations have been hampered by the difficulty of incorporat- 
ing inhomogeneities into perfectly periodic systems. Sev- 
eral theoretical approaches have been formulated to address 
this specific issue, including the virtual crystal approximation 
(VCA)^^ and coherent potential approximation (CR4).u~eI 

The VCA has been performed for over 30 years with mixed 
results; a composite potential is constructed which represents 
the average of the component atoms comprising the inhomo- 
geneity. Previous implementations of the VCA have studied 
how structural and electronic properties of semiconducting 
and ferromagnetic alloys depend on composition. These stud- 
ies typically provide qualitative and in some cases quantitative 
agreerneat with experimental and large-scale solid-state calcu- 
lations.Q|j In some^cases however, these methods have given 
unphysical results.Efu The inability of these traditional ver- 
sions of the VCA to represent differences in chemical bonding 
and ionicity accurately Jnas been suggested as the reason for 
these erroneous results. B In principle, the VCA affords great 
versatility, applicable wherever compositional inhomogeneity 
can treated in an averaged way. On its own, the VCA does not 
describe the different local atomic environment in inhomoge- 
neous materials. Local features can be studied by perform- 
ing a VCA calculation and adding in the differences between 
true atomipjjotentials and the VCA using a linear response 
method.Eiy 

For predictions of alloy densities— of states, the CPA 
has achieved considerable success.li3~li3 In this approach, a 
Green's function is constructed such that an electron encoun- 
tering either element of the alloy will not scatter on aver- 
age. Recently, nearest-neighbor charge correlations have been 
included into CR\ calculations £3 Despite this success, it is 
not convenient to perform first-principles calculations within 
the CPA which include self-consistency, Hellmann-Feynman 
forces and atomic motion. Therefore, considerable effort has 



been exerted toward the development of a VCA appropriate to 
first-principles calculations. 

In response to the shortcomings of traditional approaches 
to the VCA, we have formulated a new approach that pre- 
serves important atomic properties (atomic orbital energies, 
ionicity, polarizability, and average size of the virtual atom). 
The VCA method presented in this paper represents a consid- 
erable departure from other first-principles VCA treatments. 
A variety of calculations have been performed by averaging 
pseudopotentials of the component atoms and finding self- 
consistent wave functions and charge density in the presence 
of this average potential. As we discuss below, the nonlin- 
earity of the self-consistent Kohn-Sham equation insures that 
the average potential does not yield averaged properties. Ac- 
cordingly, we construct a potential which does give correctly 
averaged atomic eigenvalues for a wide variety of electronic 
configurations. We anticipate that our new VCA method, be- 
cause it provides correctly averaged atomic properties, will 
yield more accurate virtual atoms and provide greater agree- 
ment with superlattice solid-state calculations than traditional 
VCAs based on averaged potentials. The VCA construction 
methods are described for the combination of two atoms ac- 
cording to Ai_j;B^. In addition, we restrict our description to 
homovalent atoms. These methods can be generalized to the 
averaging of more than two atoms and heterovalency. It is im- 
portant to emphasize that since we construct a single potential 
at the atomic level, the use of our VCA pseudopotentials in 
solid-state calculations requires no additional computational 
effort. 

We have chosen to assess the accuracy of the two methods 
to reproduce superlattice calculations by computing the rel- 
ative energy differences between different phases of a feiTO- 
electric solid solution {x — 0.5) of PbZrOs and PbTiOa. This 
composition lies very near the tetragonal-rhombohedral mor- 
photropic phase boundary of this material. After establishing 
the accuracy of our new method, we use it to predict the en- 
ergetic ordering of Pb(Zri_2^Ti2:)03 (PZT) phases at another 
value of X. We have previously demonstrated the effectiveness 
of our new approach for phase stability and for strain-induced 
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phase transitions. 01n this paper, we perform a more complete 
treatment of the PZT work, involving three structural phases 
and studying a compositional phase transition. 

ADVCA (Averaged Descreened Potentials VCA): In this ap- 
proach, two independent atomic pseudopotentials (\pg and 
Vpg) are constructed. These descreened potentials are then 
averaged according to: 
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In order to express the above semilocal VCA potential in a 
nonlocal formJlH we must also know the semilocal pseudo- 
wave-functions. Therefore, we solve for the reference-state 
pseudo-wave-functions of and use these to construct 

the nonlocal potential. 

RRVCA (New Method): The ADVCA construction does not 
guarantee that the virtual atom will have atomic properties 
which are the correct average of the component atoms. There- 
fore, we impose averaged properties as a new criterion in our 
RRVCA construction. We define an averaged eigenvalue for 
the l-th state of the virtual atom. 



TABLE I. Construction parameters for the Zr, Ti and Zri-^^Tii, 
virtual crystal (VCA) designed nonlocal pseudopotentials. The Zr 
and Ti potentials were generated with the methods described in ref- 
erences |l^ and |2^. The VCA potentials were generated using the 
methods described in text. Core radii (Vc) are in atomic units, qc 
parameters are in Ry^^^, step widths are in atomic units and step 
heights are in Ry. 



Atom 






Step 
Range 


Step 
Height 


Zr(4sMp''4d") 


1.51,1.51,1.90 


7.07 


0.00-1.51 


0.77 


Ti (Ss^Sp^Sd") 


1.32,1.20,1.40 


7.07 


0.00-1.15 


5.49 


ADVCA (s^p^cf) 
X = 0.5 


generated by averaging 
Zr and Ti potentials 




RRVCA {s^p^dP) 
X = 0.5 


1.38,1.51,1.40 


7.07 


0.00-0.56 
0.56-0.79 


1.18 
1.34 


RRVCA {s^p^'d") 
X = 0.4 


1.27,1.38,1.61 


7.07 


0.00-1.25 


10.32 



^AVG 



(1 -a;)e^i +a;£^ 



(2) 



where and e^q are the valence all-electron eigenvalues for 
the nl- and n'l-th states of atoms A and B, and n need not 
equal n'. We seek a VCA potential which guarantees that 
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for the reference state. As a first step, the 
bare nuclear Coulombic potentials of the component atoms 
are averaged: 
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In addition, we determine a core charge density that is 
the weighted average of the core charge densities of the all- 
electron component atoms (/O^i-o and jO^^j.^,), 



RRVCA 
r^corc 



(1 •^)Pcorc + '''/-'core- 



(4) 



By averaging the core charge densities and all-electron po- 
tentials, we insure that the resulting virtual atom will possess 
the proper averaged size of the two component atoms. Using 
this nuclear potential and frozen-core charge density, we find 
new all-electron wave functions for the virtual atom valence 
states. We accomplish this by completing a self-consistent in- 
ward solvellj for the valence wave functions according to 
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3. J^\^f^''^^{r)\\^dr = 

where Vc is the pseudopotential core radius for the valence 
state. Operationally, we complete an inward solve well below 
Tc, insuring accurate first and second derivative determination 
for all values of r greater than Vc. These derivatives ra^e re- 
quired for the optimized pseudopotential construction.tj 



Since the form of these all-electron RRVCA wave func- 
tions will be modified within the core region when the wave 
functions are pseudizedjl3 the Kohn-Sham equations need not 
be solved between r=0 and r=rc- We construct the remain- 
der of the wave function as a smooth analytic form insuring 
norm-conservation and agreement with the inward solution. 
A new valence charge density is constructed, and the pro- 
cess is iterated to self-consistency. These wave functions are 
self-consistent solutions to the Kohn-Sham equation outside 
Tc, with eigenvalues and total norm outside rc which are ex- 
actly the average of the component atoms. Optimized semilo- 
cal pseudopotentials are then constructed for this set of wave 
functions and eigenvalues. It is important to note that any 
method for pseudopotential construction may be used because 
the averaging occurs at the all-electron level. 

The pseudopotential construction procedure guarantees the 
agreement of the RRVCA eigenvalues with the averaged all- 
electron eigenvalues for the reference state only. To im- 
prove the transferability of the RRVCA potential to other elec- 
tronic configurations, we have, employed the designed nonlo- 
cal pseudopotential method.E] This method relies on the ar- 
bitrariness of the separation between the local and semilo- 
cal coiTection terms that are used in the nonlocal potential. 
By adjusting the separation, we can dramatically improve the 
transferability of the RRVCA potential. Improving the agree- 
ment with the all-electron atomic orbital energies for a range 
of electronic configurations insures correct ionicity and po- 
larizability of the virtual atom. It is important to note that any 
method of pseudopotential construction sufficiently flexible to 
achieve eigenvalue agreement for many electronic configura- 
tions will also improve the quality of the virtual atom. 

We have appHed the ADVCA and RRVCA methods to the 
Zri-^Ti^; virtual atom at x = 0.5. All atomic energy cal- 
culations were dOjtie within the local density approximaticp, 
and the optimizedtj and designed nonlocal pseudopotentialcJ 
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TABLE II. Configuration testing for Zr, Ti and Zri-i^Tij, virtual 
crystal (VCA) atoms generated with the ADVCA and RRVCA meth- 
ods described in text. Averaged eigenvalues and total energy differ- 
ences are given for the Zr and Ti all-electron (AE) and component 
pseudopotentials (PS). Absolute errors are computed as a difference 
from the averaged component pseudopotential results. All energies 
are in Ry. 





X = 0.5 


X = 0.5 


X = 0.5 


X = 0.5 




AE 


PS 


ADVCA 


RRVCA 


State 


Energy 


Energy 


Error 


Error 




-7.5701 


-7.5701 


0.2108 


0.0000 


/ 


-5.9530 


-5.9530 


0.1554 


0.0000 




-2.5554 


-2.5581 


0.0437 


-0.0038 




-3.4488 


-3.4488 


0.2607 


0.0000 


AEtot 


0.0000 


0.0000 


0.0000 


0.0000 




-6.7832 


-6.7808 


0.1942 


0.0016 


/ 


-5.1737 


-5.1720 


0.1398 


0.0015 




-2.0292 


-2.0304 


0.0335 


-0.0009 


d° 


-2.7037 


-2.7034 


0.2362 


0.0041 


ABtot 


-2.3049 


-2.3067 


0.0387 


-0.0019 




-6.3611 


-6.3652 


0.1434 


-0.0031 


/ 


-4.7701 


-4.7744 


0.0923 


-0.0049 




-1.8390 


-1.8399 


0.0147 


0.0000 


d} 


-2.3479 


-2.3571 


0.1782 


-0.0033 




-2.8876 


-2.8936 


0.2201 


-0.0012 




-5.6569 


-5.6586 


0.1266 


-0.0029 


/ 


-4.0711 


-4.0734 


0.0764 


-0.0047 




-1.3449 


-1.3453 


0.0119 


0.001 1 




-I.68I4 


-1.6889 


0.1549 


-0.0018 


AEtot 


-4.4856 


-4.4923 


0.2348 


-0.0004 




-4.1762 


-4.1730 


0.0398 


-0.0077 


/ 


-2.6089 


-2.6064 


-0.0059 


-0.0106 




-0.3301 


-0.3293 


-0.0097 


0.0018 




-0.3205 


-0.3239 


0.0586 


-0.0072 


AStot 


-6.2574 


-6.2696 


0.3397 


-0.0037 



methods were used. The generation parameters for the com- 
ponent potentials and the VCA potentials are included in Ta- 
ble I. For all atoms, semi-core states were included as valence. 
It is important to note that although we have included multi- 
ple s-channel states, only one s nonlocal projector is included. 
For all atoms we have used the s-potential as the local po- 
tential onto which we add one or two square-step augmenta- 
tion operators. The transferability data for the ADVCA and 
RRVCA atoms are presented in Table II. For comparison we 
have also included the all-electron averaged eigenvalues and 
the averaged component pseudopotential eigenvalues. All er- 
rors in the VCA methods are computed as the difference from 
the averaged component pseudopotential eigenvalues. From 
the table it is clear that the RRVCA provides the more trans- 
ferable potentials; the RRVCA gives atomic properties which 
are orders of magnitude closer to the weighted average of the 
component atoms than the ADVCA. We now test both VCA 



potentials, using them in solid-state calculations and compar- 
ing to superlattice results. We have completed first-principles 
calculations using density functional theory and the local den- 
sity approximation on a Pb(Z]jMTio.5)03 superlattice with 
stacking in the (111) directioncil and (Oil) direction. The 
electronic wave functions are expanded in a plane-wave basis 
using a cutoff energy of 50 Ry. In addition to the semi-core 
states mentioned above, the 5d shell is included in the Pb po- 
tential. Scalar relativistic effects are included in the generation 
of the Pb pseudopotential. Brillouin zone integrations are ap- 
proximated accurately as sums on a 4x4x4 Monkhorst-Pack 
/c-point meshB3 

The experimentally determined phase diagram for PZT 
shows two phase boundaries relevant to the present work. 
A nearly temperature-independent boundary at x ^ 0.45 
separates the Ti-rich tetragonal PAmm phase, from the Zr- 
rich rhombohedral (i?3c and Rira) phases. Around 400- 
500K, the rhombohedral region exhibits a boundary between 
i?3c (low-temperature) and R'dm (high-leuuperature) phases, 
which depends weakly on compositionBScj The i?3c phase 
shows complex oxygen octahedral tilting, which doubles the 
primitive unit cell to ten atoms. EZl 

Table III contains the x — 0.5 PZT superlattice equilib- 
rium cell dimensions and energy differences for the tetrago- 
nal, low- and high-temperature rhombohedral phases for both 
cation orderings. For all rhombohedral calculations, we have 
neglected the small shear relaxations. Energy differences have 
been computed relative to the tetragonal phase. For the (111) 
ordering, we find that the low- and high-temperature rhombo- 
hedral phases are 0.035 eV and 0.05 eV higher in energy than 
the tetragonal phase. For the (01 1) ordering, we find the low- 
and high-temperature rhombohedral phases are 0.093 eV and 
0.115 eV above the tetragonal ground state. For both order- 
ings, the low-temperature rhombohedral phase is lower than 
the high-temperature phase. 

For both VCA methods, we have completed full elec- 
tronic and structural relaxations for five-atom unit cells which 
possess tetragonal and _R3to symmetries. For the low- 
temperature rhombohedral phase, we have studied a 10-atom 
unit cell which possesses _R3c symmetry. In Table III, we 
report the results of the ADVCA. ADVCA predicts that the 
ground-state phase for this composition of PZT is the high- 
temperature rhombohedral phase, in direct opposition to the 
superlattice calculations and experimental observations. In 
addition we find that ADVCA significantly underestimates the 
c/a ratio for the tetragonal phase. We attribute the lack of 
agreement to the poor quality of the ADVCA virtual atom. 
The second to last column of Table III contains the results 
from the RRVCA. We find the same energy ordering as the 
superlattice calculations and experiment, with the low- and 
high-temperature phases 0.05 eV and 0.072 eV higher than 
the tetragonal phase. Finally, we find excellent agreement in 
the structural parameters between the RRVCA and superlat- 
tice calculations. We attribute the quantitative agreement in 
the energetic and structural results to the improved quality of 
the RRVCA virtual atom. 

Having demonstrated the effectiveness of the RRVCA and 
the inability of the ADVCA to reproduce the superlattice re- 
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TABLE III. Solid-state calculation results for superlattice and 
virtual crystals of x = 0.5 and x = 0.4 compositions of 
Pb(Zri_iTia:)03 with (111) and (Oil) orderings. Energy differences 
are relative to the tetragonal phase in all cases. Energies and struc- 
tural parameters are given for a 40-atom unit cell. 

Superlattice 
x^O.5 x = 0.5 x = 0.5 x = 0.5 x = 0.4 
Structure (111) (Oil) ADVCA RRVCA RRVCA 



Rhombohedral 


a 


8.052 


8.029 


8.077 


8.043 


8.033 


{mm) 


AE 


0.050 


0.115 


-0.041 


0.072 


0.048 


Tetragonal 


c 


8.199 


8.219 


8.131 


8.210 


8.199 


{PAmm) 


a 


7.869 


7.989 


8.054 


7.952 


7.957 




c/a 


1.041 


1.029 


1.010 


1.033 


1.030 




AE 


0.000 


0.000 


0.000 


0.000 


0.000 


Rhombohedral 


a 


8.026 


8.004 


8.046 


8.007 


8.019 


(i?3c) 


AE 


0.035 


0.093 


-0.067 


0.050 


-0.044 


suits for X = 


0.5, 


we use 


the RRVCA to 


examine 


further 



the compositional phase transition in PZT. We have chosen 
to study the x = 0.4 composition of PZT, since the ground- 
state structure for this composition lies in the low-temperature 
rhombohedral phase region of the phase diagram. We have 
constructed a x ^ 0.4 Zri_j^Tij^ potential using the RRVCA 
according to the parameters in Table I. The transferability of 
this potential is equal to that of the x = 0.5 RRVCA po- 
tential. We have completed structural relaxations of the two 
five-atom unit cell structures described above as well as the 
ten-atom low-temperature rhombohedral phase. The last col- 
umn of Table III contains the structural parameters for the 
three phases and energy differences relative to the tetrago- 
nal phase. We find that moving from x = 0.5 to x — 0.4 
composition shifts both rhombohedral phases lower in energy 
relative to the tetragonal phase. The high-temperature rhom- 
bohedral phase energy drops by approximately 0.012 eV and 
the low-temperature phase decreases by 0.094 eV, becoming 
the ground-state structure for the x — 0.4 composition. This 
result is in direct agreement with the experimental findings 
and demonstrates the ability of the RRVCA to predict compo- 
sitional phase transitions. 

In this paper we have presented a new method for con- 
structing virtual crystal pseudopotentials and have applied this 
method to the Ti and Zr atoms. Our new method is based 
on the self-consistent determination of averaged all-electron 
properties which insures exact agreement of the proper aver- 
aged all-electron eigenvalues for the reference configuration 
of the virtual atom. This procedure also yields the proper av- 
eraged atomic size. To improve the electronic properties of 
the virtual atom at all other electronic configurations, we em- 
ployed the designed nonlocal pseudopotential approach. We 
have compared our new method to a more traditional one in 
which the descreened pseudopotentials are averaged, comput- 
ing the relative energy differences for the high-temperature 
rhombohedral, low-temperature rhombohedral and tetragonal 
phases of Pb(Zri_a;Tia;)03 and comparing those results to su- 



perlattice calculations with two different cationic orderings. 
For this more traditional method, we find neither quantitative 
nor qualitative agreement with superlattice results. However, 
for our new method, we find excellent agreement not only in 
the energetics of the three phases but also the structural pa- 
rameters. Finally we use our new method to predict the com- 
positional phase transition between x = 0.5 and x = 0.4 
Pb(Zri_2;Ti2;)03. Our new method predicts that the low- 
temperature rhombohedral phase is the ground-state structure 
for the X = 0.4 composition in agreement with experiment. 
This calculation represents the first ab initio determination of 
a compositional phase boundary in an ferroelectric oxide. 

We have recently become aware of an alternative approach 
to imncDving the first-principles VCA. Bellaiche and Van- 
derbilt£3 apply the nonlocal projectors of all the component 
atoms at the_site of the virtual atom. Using ultrasoft pseu- 
dopotentials,Ej they find accurate piezoelectric and ferroelec- 
tric properties for Pb(Zro.5Tio.5)03 using their approach. 
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